The goal of this analysis is to determine the enrichment of m6A-Tracer signal around nucleoli (marked by Ki67), following Ki67 pA-DamID. Dam-only signal is the control.
Load the required libraries.
# Load libraries
library(ggplot2)
library(tidyverse)
library(ggbeeswarm)
library(pixmap)
# Input data
results_dir <- c("../ts201008_E1304_E1321_confocal_various/ts201008_pADamID_m6ATracer_analysis/",
"../ts201102_E1349_confocal_synchronization/ts201102_RPE_m6ATracer/ts201102_RPE_timecourse_analysis/",
"../ts210308_E1499_confocal/ts210316_RPE_m6ATracer/ts210308_RPE_m6ATracer_analysis/")
####"../ts210402_E1546_confocal_Ki67_antibodies/ts210402_HCT116_wt_pADamID_Ki67_antibodies/ts210402_HCT116_wt_pADamID_Ki67_antibodies_analysis/")
# Prepare output directory
output_dir <- "ts220503_nucleolar_enrichment_cell_types"
dir.create(output_dir, showWarnings = FALSE)
# Prepare output fileslibrary(knitr)
opts_chunk$set(dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/"),
cache = T)
pdf.options(useDingbats = FALSE)ReadPGN <- function(file.name) {
# Get image, only grey values
suppressWarnings(read.pnm(file.name)@grey*255)
}The main problem is substantial difference in microscopy slide quality between the replicates / cell types. This results in consistent results, but also in qualitative variation. I will explore a bit whether I can solve this, but otherwise I can always just “explain” this in the text. Solving means tweaking the parameters to segment nucleoli, mostly.
List the different samples.
metadata <- tibble(cell_statistics = dir(results_dir, recursive = T,
full.names = T,
pattern = "cells_statistics")) %>%
mutate(sample = str_remove(basename(cell_statistics),
"_cells_.*"),
sample = str_remove(sample, "r._"),
sample = str_replace(sample, "unsync", "RPE")) %>%
filter(! grepl("Act", sample) &
! grepl("Preview", sample) &
! grepl("t\\d+h+_", sample) &
! grepl("NOV|HPA", sample)) %>%
mutate(replicate = case_when(grepl("r1_", cell_statistics) ~ "r2",
grepl("r2_", cell_statistics) ~ "r3",
T ~ "r1"),
replicate = factor(replicate, levels = c("r1", "r2", "r3")),
cell = case_when(grepl("HCT116", sample) ~ "HCT116",
grepl("K562", sample) ~ "K562",
T ~ "RPE"),
cell = factor(cell, levels = c("RPE", "HCT116", "K562")),
condition = case_when(grepl("DMSO", cell_statistics) ~ "DMSO",
T ~ "wt"),
target = ifelse(grepl("Ki67", sample), "Ki67", "Dam"),
target = factor(target, levels = c("Dam", "Ki67")),
slide = str_remove(sample, ".*_")) %>%
mutate(sample = paste(cell, target, condition, replicate, slide, sep = "_")) %>%
mutate(objects = str_replace(cell_statistics,
"cells_statistics.csv",
"dapi_segment.pgm"),
cell_internal = str_replace(objects,
"dapi_segment",
"dapi_internal"),
nucleolus_internal = str_replace(objects,
"dapi_segment",
"nucleolus_internal"),
nucleolus_external = str_replace(objects,
"dapi_segment",
"nucleolus_external"),
nucleolus = str_replace(objects,
"dapi_segment",
"nucleolus_smooth"),
tracer = str_replace(objects,
"dapi_segment",
"m6ATracer_smooth")) %>%
arrange(cell, target, condition, replicate, slide)
# Print metadata
metadata %>%
print(n = Inf)## # A tibble: 73 × 13
## cell_statistics sample replicate cell condition target slide objects
## <chr> <chr> <fct> <fct> <chr> <fct> <chr> <chr>
## 1 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 002 ../ts2011…
## 2 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 004 ../ts2011…
## 3 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 006 ../ts2011…
## 4 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 009 ../ts2011…
## 5 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 011 ../ts2011…
## 6 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 013 ../ts2011…
## 7 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 015 ../ts2011…
## 8 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 017 ../ts2011…
## 9 ../ts201102_E1349_… RPE_Da… r1 RPE wt Dam 019 ../ts2011…
## 10 ../ts210308_E1499_… RPE_Da… r2 RPE wt Dam 004 ../ts2103…
## 11 ../ts210308_E1499_… RPE_Da… r2 RPE wt Dam 007 ../ts2103…
## 12 ../ts210308_E1499_… RPE_Da… r2 RPE wt Dam 009 ../ts2103…
## 13 ../ts210308_E1499_… RPE_Da… r2 RPE wt Dam 012 ../ts2103…
## 14 ../ts210308_E1499_… RPE_Da… r2 RPE wt Dam 015 ../ts2103…
## 15 ../ts210308_E1499_… RPE_Da… r3 RPE wt Dam 002 ../ts2103…
## 16 ../ts210308_E1499_… RPE_Da… r3 RPE wt Dam 005 ../ts2103…
## 17 ../ts210308_E1499_… RPE_Da… r3 RPE wt Dam 007 ../ts2103…
## 18 ../ts210308_E1499_… RPE_Da… r3 RPE wt Dam 009 ../ts2103…
## 19 ../ts210308_E1499_… RPE_Da… r3 RPE wt Dam 011 ../ts2103…
## 20 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 009 ../ts2011…
## 21 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 011 ../ts2011…
## 22 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 014 ../ts2011…
## 23 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 016 ../ts2011…
## 24 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 018 ../ts2011…
## 25 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 021 ../ts2011…
## 26 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 023 ../ts2011…
## 27 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 025 ../ts2011…
## 28 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 027 ../ts2011…
## 29 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 029 ../ts2011…
## 30 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 033 ../ts2011…
## 31 ../ts201102_E1349_… RPE_Ki… r1 RPE wt Ki67 035 ../ts2011…
## 32 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 004 ../ts2103…
## 33 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 006 ../ts2103…
## 34 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 008 ../ts2103…
## 35 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 010 ../ts2103…
## 36 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 013 ../ts2103…
## 37 ../ts210308_E1499_… RPE_Ki… r2 RPE wt Ki67 015 ../ts2103…
## 38 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 002 ../ts2103…
## 39 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 004 ../ts2103…
## 40 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 006 ../ts2103…
## 41 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 008 ../ts2103…
## 42 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 010 ../ts2103…
## 43 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 012 ../ts2103…
## 44 ../ts210308_E1499_… RPE_Ki… r3 RPE wt Ki67 015 ../ts2103…
## 45 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Dam 003 ../ts2010…
## 46 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Dam 005 ../ts2010…
## 47 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Dam 007 ../ts2010…
## 48 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Dam 009 ../ts2010…
## 49 ../ts201008_E1304_… HCT116… r1 HCT1… wt Dam 003 ../ts2010…
## 50 ../ts201008_E1304_… HCT116… r1 HCT1… wt Dam 006 ../ts2010…
## 51 ../ts201008_E1304_… HCT116… r1 HCT1… wt Dam 008 ../ts2010…
## 52 ../ts201008_E1304_… HCT116… r1 HCT1… wt Dam 010 ../ts2010…
## 53 ../ts201008_E1304_… HCT116… r1 HCT1… wt Dam 012 ../ts2010…
## 54 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 003 ../ts2010…
## 55 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 005 ../ts2010…
## 56 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 007 ../ts2010…
## 57 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 008 ../ts2010…
## 58 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 009 ../ts2010…
## 59 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 011 ../ts2010…
## 60 ../ts201008_E1304_… HCT116… r1 HCT1… DMSO Ki67 012 ../ts2010…
## 61 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 005 ../ts2010…
## 62 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 007 ../ts2010…
## 63 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 009 ../ts2010…
## 64 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 011 ../ts2010…
## 65 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 013 ../ts2010…
## 66 ../ts201008_E1304_… K562_D… r1 K562 wt Dam 015 ../ts2010…
## 67 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 004 ../ts2010…
## 68 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 006 ../ts2010…
## 69 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 008 ../ts2010…
## 70 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 010 ../ts2010…
## 71 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 014 ../ts2010…
## 72 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 016 ../ts2010…
## 73 ../ts201008_E1304_… K562_K… r1 K562 wt Ki67 018 ../ts2010…
## # … with 5 more variables: cell_internal <chr>, nucleolus_internal <chr>,
## # nucleolus_external <chr>, nucleolus <chr>, tracer <chr>
Load the cells and show the sizes. Possibly filter on those.
# Load all cells
LoadCells <- function(i) {
tib <- read_csv(metadata$cell_statistics[i]) %>%
dplyr::select(c(1:2)) %>%
rename_all(~c("volume", "surface")) %>%
add_column(object = 1:nrow(.),
sample = metadata$sample[i],
replicate = metadata$replicate[i],
cell = metadata$cell[i],
condition = metadata$condition[i],
target = metadata$target[i],
slide = metadata$slide[i]) %>%
mutate(cell_id = paste(sample, object, sep = "_"))
tib
}
tib_cells <- lapply(1:nrow(metadata), LoadCells) %>%
purrr::reduce(bind_rows)Plot and filter these cells.
# Set size limits
limits_size <- c(10, 125)
# Plot cell volume
tib_cells %>%
ggplot(aes(x = target, y = volume)) +
geom_hline(yintercept = 0, alpha = 0) +
geom_hline(yintercept = limits_size, col = "black", linetype = "dashed") +
geom_quasirandom(col = "darkgrey") +
geom_boxplot(fill = NA, col = "black", outlier.shape = NA) +
facet_grid(. ~ cell) +
xlab("") +
ylab("Volume slice (micron)") +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# Filter cells for cutoff
tib_cells <- tib_cells %>%
mutate(size_filter = volume >= limits_size[1] &
volume <= limits_size[2]) Load the .pgm files for the various channels.
# Load all cells
LoadPGMMaps <- function(i, n_min = 5) {
# Load objects
objects <- ReadPGN(metadata$objects[i])
nucleolus <- ReadPGN(metadata$nucleolus[i])
nucleolus_internal <- ReadPGN(metadata$nucleolus_internal[i])
nucleolus_external <- ReadPGN(metadata$nucleolus_external[i])
tracer <- ReadPGN(metadata$tracer[i])
tib <- tibble(object = c(objects),
nucleolus_internal = c(nucleolus_internal),
nucleolus_external = c(nucleolus_external),
nucleolus = c(nucleolus),
tracer = c(tracer)) %>%
filter(object != 0) %>%
rowwise() %>%
mutate(nucleolus_internal = min(nucleolus_internal, 254),
nucleolus_external = min(nucleolus_external, 254)) %>%
ungroup() %>%
group_by(object, nucleolus_internal, nucleolus_external) %>%
summarise(n = n(),
tracer = mean(tracer),
nucleolus = mean(nucleolus)) %>%
ungroup() %>%
group_by(object) %>%
mutate(tracer_enrichment = log2(tracer / weighted.mean(tracer, n)),
nucleolus_enrichment = log2(nucleolus / weighted.mean(nucleolus, n))) %>%
ungroup() %>%
filter(n > n_min) %>%
mutate(distance = nucleolus_external - nucleolus_internal) %>%
arrange(object, distance) %>%
add_column(sample = metadata$sample[i],
replicate = metadata$replicate[i],
cell = metadata$cell[i],
condition = metadata$condition[i],
target = metadata$target[i],
slide = metadata$slide[i]) %>%
mutate(cell_id = paste(sample, object, sep = "_"))
tib
}
tib_intensities <- lapply(1:nrow(metadata), LoadPGMMaps) %>%
purrr::reduce(bind_rows)
# Filter out cells that did not pass the cell size filter
tib_intensities <- tib_intensities %>%
dplyr::filter(cell_id %in% (tib_cells %>%
filter(size_filter == T) %>%
pull(cell_id)))Plots.
# Filter very bad cells
# - only positive / negative scores
# - extreme values
remove_samples <- tib_intensities %>%
gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
mutate(key = str_remove(key, "_enrichment")) %>%
filter(distance > -11 & distance <= 11) %>%
group_by(key, target, cell, cell_id) %>%
dplyr::summarise(q1 = quantile(value, 0.1),
q5 = quantile(value, 0.5),
q9 = quantile(value, 0.9)) %>%
filter(q1 > 0.5 | q9 < 0 | q1 < -2.5 | q9 > 3) %>%
pull(cell_id)## `summarise()` has grouped output by 'key', 'target', 'cell'. You can override using the `.groups` argument.
# Plot intensity per cell
tib_gather <- tib_intensities %>%
filter(! cell_id %in% remove_samples) %>%
gather(key, value, nucleolus, tracer) %>%
filter(distance > -14 & distance <= 14) %>%
mutate(key = factor(key, levels = c("tracer", "nucleolus")))
tib_gather %>%
ggplot(aes(x = distance, y = value, col = key,
group = interaction(key, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line() +
facet_grid(target ~ cell, scales = "free_y") +
xlab("Distance to nucleolus (pixel)") +
ylab("Intensity (A.U.)") +
coord_cartesian(xlim = c(-10, 10)) +
scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# Plot intensity per cell - for enrichments
tib_gather <- tib_intensities %>%
filter(! cell_id %in% remove_samples) %>%
gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
filter(distance > -14 & distance <= 14) %>%
mutate(key = str_remove(key, "_enrichment")) %>%
mutate(key = factor(key, levels = c("tracer", "nucleolus")))
tib_gather %>%
ggplot(aes(x = distance, y = value, col = key,
fill = target, group = interaction(key, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line(alpha = 0.1) +
stat_summary(fun = mean, aes(group = key),
geom = "line", size = 2) +
facet_grid(cell ~ target, scales = "free_y") +
xlab("Distance to nucleolus (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(-10, 10),
ylim = c(-2.1, 2.3)) +
scale_color_manual(values = c("grey30", "red")) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))tib_gather %>%
ggplot(aes(x = distance, y = value, col = target,
fill = target,
group = interaction(target, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line(alpha = 0.1) +
stat_summary(fun = mean, aes(group = target),
geom = "line", size = 2) +
# stat_summary(fun.data = mean_se, geom = "ribbon", aes(group = target),
# fun.args = list(mult = 1.96), col = NA, alpha = 0.2) +
facet_grid(cell ~ key, scales = "free_y") +
xlab("Distance to nucleolus (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(-10, 10),
ylim = c(-2.1, 2.3)) +
scale_color_manual(values = c("grey30", "red")) +
scale_fill_manual(values = c("grey30", "red")) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# Means only
# 1) Combined
tib_gather %>%
ggplot(aes(x = distance, y = value, col = target,
fill = target, group = target)) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
stat_summary(fun = mean, geom = "line", size = 1) +
# stat_summary(fun.data = mean_se, geom = "errorbar",
# col = "darkgrey", linetype = "solid") +
stat_summary(fun.data = mean_se, geom = "ribbon",
col = NA, alpha = 0.2) +
facet_grid(cell ~ key) +
xlab("Distance to nucleolus (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(-10, 10),
ylim = c(-1.5, 1.9)) +
#scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# 2) Per replicate
tib_gather %>%
ggplot(aes(x = distance, y = value, col = target, linetype = interaction(replicate, condition),
fill = target, group = interaction(target, replicate, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
stat_summary(fun = mean, geom = "line", size = 1) +
# stat_summary(fun.data = mean_se, geom = "errorbar",
# col = "darkgrey", linetype = "solid") +
stat_summary(fun.data = mean_se, geom = "ribbon",
col = NA, alpha = 0.2) +
facet_grid(cell ~ key) +
xlab("Distance to nucleolus (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(-10, 10),
ylim = c(-1.5, 1.9)) +
#scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))As requested by reviewer #2, include quantification from periphery of nuclei
# Load all cells
LoadPGMMaps <- function(i, n_min = 5) {
# Load objects
objects <- ReadPGN(metadata$objects[i])
cell_internal <- ReadPGN(metadata$cell_internal[i])
nucleolus <- ReadPGN(metadata$nucleolus[i])
tracer <- ReadPGN(metadata$tracer[i])
tib <- tibble(object = c(objects),
cell_internal = c(cell_internal),
nucleolus = c(nucleolus),
tracer = c(tracer)) %>%
filter(object != 0) %>%
rowwise() %>%
mutate(cell_internal = min(cell_internal, 254)) %>%
ungroup() %>%
group_by(object, cell_internal) %>%
summarise(n = n(),
tracer = mean(tracer),
nucleolus = mean(nucleolus)) %>%
ungroup() %>%
group_by(object) %>%
mutate(tracer_enrichment = log2(tracer / weighted.mean(tracer, n)),
nucleolus_enrichment = log2(nucleolus / weighted.mean(nucleolus, n))) %>%
ungroup() %>%
filter(n > n_min) %>%
mutate(distance = 254 - cell_internal) %>%
arrange(object, distance) %>%
add_column(sample = metadata$sample[i],
replicate = metadata$replicate[i],
cell = metadata$cell[i],
condition = metadata$condition[i],
target = metadata$target[i],
slide = metadata$slide[i]) %>%
mutate(cell_id = paste(sample, object, sep = "_"))
tib
}
tib_intensities <- lapply(1:nrow(metadata), LoadPGMMaps) %>%
purrr::reduce(bind_rows)
# Filter out cells that did not pass the cell size filter
tib_intensities <- tib_intensities %>%
dplyr::filter(cell_id %in% (tib_cells %>%
filter(size_filter == T) %>%
pull(cell_id)))Plots.
# Filter same bad cells as before
# Plot intensity per cell
tib_gather <- tib_intensities %>%
filter(! cell_id %in% remove_samples) %>%
gather(key, value, nucleolus, tracer) %>%
filter(distance <= 14) %>%
mutate(key = factor(key, levels = c("tracer", "nucleolus")))
tib_gather %>%
ggplot(aes(x = distance, y = value, col = key,
group = interaction(key, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line() +
facet_grid(target ~ cell, scales = "free_y") +
xlab("Distance from periphery (pixel)") +
ylab("Intensity (A.U.)") +
coord_cartesian(xlim = c(0, 10)) +
scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# Plot intensity per cell - for enrichments
tib_gather <- tib_intensities %>%
filter(! cell_id %in% remove_samples) %>%
gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
filter(distance > -14 & distance <= 14) %>%
mutate(key = str_remove(key, "_enrichment")) %>%
mutate(key = factor(key, levels = c("tracer", "nucleolus")))
tib_gather %>%
ggplot(aes(x = distance, y = value, col = key,
fill = target, group = interaction(key, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line(alpha = 0.1) +
stat_summary(fun = mean, aes(group = key),
geom = "line", size = 2) +
facet_grid(cell ~ target, scales = "free_y") +
xlab("Distance from periphery (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(0, 10),
ylim = c(-3, 0.5)) +
scale_color_manual(values = c("grey30", "red")) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))tib_gather %>%
ggplot(aes(x = distance, y = value, col = target,
fill = target,
group = interaction(target, cell_id, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
geom_line(alpha = 0.1) +
stat_summary(fun = mean, aes(group = target),
geom = "line", size = 2) +
# stat_summary(fun.data = mean_se, geom = "ribbon", aes(group = target),
# fun.args = list(mult = 1.96), col = NA, alpha = 0.2) +
facet_grid(cell ~ key, scales = "free_y") +
xlab("Distance from periphery (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(0, 10),
ylim = c(-3, 0.5)) +
scale_color_manual(values = c("grey30", "red")) +
scale_fill_manual(values = c("grey30", "red")) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# Means only
# 1) Combined
tib_gather %>%
ggplot(aes(x = distance, y = value, col = target,
fill = target, group = target)) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
stat_summary(fun = mean, geom = "line", size = 1) +
# stat_summary(fun.data = mean_se, geom = "errorbar",
# col = "darkgrey", linetype = "solid") +
stat_summary(fun.data = mean_se, geom = "ribbon",
col = NA, alpha = 0.2) +
facet_grid(cell ~ key) +
xlab("Distance from periphery (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(0, 10),
ylim = c(-2, 0.2)) +
#scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))# 2) Per replicate
tib_gather %>%
ggplot(aes(x = distance, y = value, col = target, linetype = interaction(replicate, condition),
fill = target, group = interaction(target, replicate, condition))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
stat_summary(fun = mean, geom = "line", size = 1) +
# stat_summary(fun.data = mean_se, geom = "errorbar",
# col = "darkgrey", linetype = "solid") +
stat_summary(fun.data = mean_se, geom = "ribbon",
col = NA, alpha = 0.2) +
facet_grid(cell ~ key) +
xlab("Distance from periphery (pixel)") +
ylab("Enrichment over mean (log2)") +
coord_cartesian(xlim = c(0, 10),
ylim = c(-2, 0.2)) +
#scale_color_grey() +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))Several observations:
Overall, this is a positive result.
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.36 pixmap_0.4-12 ggbeeswarm_0.6.0 forcats_0.5.1
## [5] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.1.0
## [9] tidyr_1.1.4 tibble_3.1.6 tidyverse_1.3.1 ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 lubridate_1.8.0 assertthat_0.2.1 digest_0.6.28
## [5] utf8_1.2.2 R6_2.5.1 cellranger_1.1.0 backports_1.4.0
## [9] reprex_2.0.1 evaluate_0.14 highr_0.9 httr_1.4.2
## [13] pillar_1.6.4 rlang_0.4.12 readxl_1.3.1 rstudioapi_0.13
## [17] jquerylib_0.1.4 rmarkdown_2.11 labeling_0.4.2 bit_4.0.4
## [21] munsell_0.5.0 broom_0.7.10 compiler_4.1.2 vipor_0.4.5
## [25] modelr_0.1.8 xfun_0.28 pkgconfig_2.0.3 htmltools_0.5.2
## [29] tidyselect_1.1.1 codetools_0.2-18 fansi_0.5.0 crayon_1.4.2
## [33] tzdb_0.2.0 dbplyr_2.1.1 withr_2.4.2 grid_4.1.2
## [37] jsonlite_1.7.2 gtable_0.3.0 lifecycle_1.0.1 DBI_1.1.1
## [41] magrittr_2.0.1 scales_1.1.1 vroom_1.5.6 cli_3.1.0
## [45] stringi_1.7.5 farver_2.1.0 fs_1.5.0 xml2_1.3.2
## [49] ellipsis_0.3.2 generics_0.1.1 vctrs_0.3.8 tools_4.1.2
## [53] bit64_4.0.5 glue_1.5.0 beeswarm_0.4.0 hms_1.1.1
## [57] parallel_4.1.2 fastmap_1.1.0 yaml_2.2.1 colorspace_2.0-2
## [61] rvest_1.0.2 haven_2.4.3